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1 ABSTRACT 

LSPRAY-II is a Lagrangian spray solver devel- 
oped for application with parallel computing and un- 
structured grids. It is designed to be massively paral- 
lel and could easily be coupled with any existing gas- 
phase flow and/or Monte Carlo Probability Density 
Function (PDF) solvers. The solver accommodates 
the use of an unstructured mesh with mixed elements 
of either triangular, quadrilateral, and/or tetrahedral 
type for the gas flow grid representation. It is mainly 
designed to predict the flow, thermal and transport 
properties of a rapidly vaporizing spray because of 
its importance in aerospace application. The manual 
provides the user with an understanding of various 
models involved in the spray formulation, its code 
structure and solution algorithm, and various other 
issues related to parallelization and its coupling with 
other solvers. With the development of LSPRAY-II, 
we have advanced the state-of-the-art in spray com- 
putations in several important ways. Some of the 
major features of the spray module are: 

• It facilitates the use of both unstructured grids 
and parallel computing and, thereby, facilitates 
large-scale combustor computations involving 
complex geometrical configurations. 

• In order to deal with modern gas-turbine fuels 
that are mixtures of many compounds, we have 
extended the spray formulation to the modeling 
of multicomponent liquid fuels. 

• In order to extend the applicability of the spray 
computations over a wide range of low-pressure 
conditions, we have completed the implemen- 
tation of the liquid-phase variable thermo- 
transport properties into our spray formulation. 

• The spray module is designed in such a way so 
that it could easily be coupled with other CFD 
codes. 


• It can be used in the calculation of both steady 
as well as unsteady computations. 

• We have developed and implemented a time- 
averaging method into the calculation of the 
liquid-phase source contributions in order to ac- 
celerate convergence towards a steady state so- 
lution. 

• The spray module has a multi-liquid and multi- 
injection capability. 

With the aim of improving the overall solution 
procedure involving sprays, we have made several rel- 
evant contributions to the gas side of the computa- 
tions: 

• In order to demonstrate the importance of chem- 
istry/turbulence interactions in the modeling 
of reacting sprays, we have extended the joint 
scalar Monte Carlo PDF (Probability Density 
Function) approach to the modeling of spray 
flames. 

• In order to account for the nonideal gas behav- 
ior under critical and supercritical conditions, 
we have completed the integration of a high- 
pressure EOS (Equation-of-State) into our CFD 
module and, also, the implementation of a high- 
pressure correction into the calculation of the 
gas-phase transport properties. 

This modeling approach provided favorable re- 
sults when applied to several different spray flames 
representative of those encountered in both gas- 
turbine combustors as well as stratified-charge rotary 
combustion (Wankel) engines. 
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2 NOMENCLATURE 


a n outward area normal vector 

of the nth surface, m 2 
Bk Spalding transfer number 

Cd drag coefficient 

C p specific heat, J/(kg K) 

D turbulent diffusion coefficient, m 2 /s 

d droplet diameter, m 

h specific enthalpy, J/kg 

k turbulence kinetic energy, m 2 /s 2 

kij binary interaction parameter 

Ik mixture latent heat of evaporation, J/kg 

lk, e ff effective latent heat of evaporation, 

J /kg (defined in Eq. (35)) 
l kin heat of vaporization at normal 

boiling point, J/kg 
Mi molecular weight of ith 

species, kg/kg-mole 
rrik droplet vaporization rate, kg/s 

rriko initial mass flow rate associated 

with kth droplet group, kg/s 
Nf number of surfaces contained in 

a given computational cell 
N p total number of computational cells 
rik number of droplets in kth group 

P pressure, N/m 2 

P c critical pressure, N/m 2 

P r Prandtl number 

R u gas constant, J/(kg K) 

Re Reynolds number 

rk droplet radius, m 

rko initial droplet radial location, m 
Sk droplet radius- squared ( = r^), m 2 
Smic liquid source contribution of the 
gas-phase continuity equation 
Smie liquid source contribution of the 
gas-phase energy equation 
Smim liquid source contribution of the 
gas-phase momentum equations 
s m i s liquid source contribution of the 
gas-phase species equations 
T temperature, K 

T n b normal boiling point, K 

T c critical temperature, K 

t time, s 

Ui ith velocity component, m/s 

Uik ith velocity component of 

kth droplet group, m/s 
V c volume of the computational cell, m 3 
or critical molar volume, m 3 /kg-mole 
V n molar volume at normal pressure, 
m 3 /kg-mole 


Wj gas-phase chemical reaction rate, 1/s 
Xi Cartesian coordinate in the ith direction, m 
yj mass fraction of jth species 

x spatial vector 

a represents a coordinate related to a Hill’s 
Vortex or a spray cone angle, deg. 

/? a spray cone angle, deg. 

X mole fraction 

A tf local time step used in the flow solver, s 

A t g i global time step used in the spray solver, s 

A tu fuel injection time step, s 

A t m i allowable time step in the spray solver, s 

AV computational cell volume, m 3 

S Dirac-delta function 

e rate of turbulence dissipation, m 2 /s 3 

€j fractional mass evaporating rate of species 

at the droplet surface 

r* turbulent diffusion coefficient, kg/ms 

A thermal conductivity, J/(ms K) 

fj, dynamic viscosity, kg/ms 

cj turbulence frequency, 1/s 

or acentric factor in Table 1 

p density, kg/m 3 

pin liquid density (at 1 bar, 273.15 K), kg/m 3 
a characteristic diameter of a molecule in Table 1 
r stress tensor term, kg/ms 2 

9 void fraction 

or a spray cone angle, deg. 

Subscripts 

/ represents conditions associated with fuel 
g global or gas-phase 

i index for the coordinate or species components 
j index for the species component 
k droplet group or liquid-phase conditions 
l liquid-phase 

n nth-face of the computational cell 
o initial conditions or oxidizer 
p conditions associated with the 
properties of a grid cell 
s represents conditions at the droplet 
surface or adjacent computational cell 
t conditions associated with time 

, partial differentiation with respect 

to the variable followed by it 

Superscripts 

n fluctuations 
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3 INTRODUCTION 

There are many occurrences of sprays in a va- 
riety of industrial and power applications and mate- 
rials processing [1]. A liquid spray is a two phase 
flow with the gas as the continuous phase and the 
liquid as the dispersed phase in the form of droplets 
or ligaments [1]. The interaction between the two 
phases, which are coupled through exchanges of mass, 
momentum, and energy, can occur in different ways 
at disparate time and length scales involving various 
thermal, mass, and fluid dynamic factors. 

A number of finite-difference formulations have 
been advanced over the years for predicting the flow 
(mass and momentum) and thermal properties of 
a rapidly vaporizing spray. Some of the pros and 
cons of various formulations can be found in [1- 
5]. Depending on the nature of the spray, an ap- 
propriate selection could be made from the choice 
of various well-known spray formulations (multicon- 
tinua, discrete-particle, or probabilistic) based on ei- 
ther a Lagrangian or an Eulerian representation for 
the liquid-phase equations by incorporating appropri- 
ately chosen droplet sub-grid models. In this manual, 
we only summarize the salient aspects of the spray 
formulation adopted from our previous work [6-15] 
without attempting to provide an in-depth review on 
the transport and fluid dynamic behavior of react- 
ing sprays. The present solution procedure could be 
used within the context of both the multicontinua and 
probabilistic spray formulations, and it allows for res- 
olution on a scale greater than the average spacing 
between two neighboring droplets [1]. An Eulerian 
scheme is assumed for the gas phase equations as it 
is an adopted choice for NCC (National Combustion 
Code). The liquid-phase equations form a system 
of hyperbolic equations and they could be solved by 
means of either an Eulerian or a Lagrangian represen- 
tation. A Lagrangian scheme is chosen as it reduces 
the errors associated with numerical diffusion. The 
droplet sub-grid models axe based on well established 
models for droplet drag; the vaporization models of a 
poly disperse spray take into account the transient ef- 
fects associated with the droplet internal heating and 
the forced convection effects associated with droplet 
internal circulation; and it employs models for gas- 
film valid over a wide range of low to intermediate 
droplet Reynolds numbers [7]. The present formula- 
tion is based on a stochastic particle tracking method 
applicable for flows with a dilute spray approxima- 
tion where the droplet loading is low. The numeri- 


cal method could be used within the context of both 
steady and unsteady calculations [8-13]. Not consid- 
ered in the present release of the code are the effects 
associated with the droplet breakup, droplet/shock 
interaction, and dense spray effects. 

In order to facilitate large-scale combustor com- 
putations, we have extended our previous work on 
sprays to parallel computing [10-14]. But it is also 
well known that considerable effort usually goes into 
generating computational grids of practical combus- 
tor geometries. In order to allow representation of 
complex geometries with relative ease, we have ex- 
tended our previous work on sprays to unstructured 
meshes [11-14] as it is well known that the grid gen- 
eration time could be improved considerably with 
the use of widely available grid generation packages 
such as GRIDGEN. With the aim of advancing the 
current multi-dimensional computational tools used 
in the design of advanced technology combustors, 
two new computer codes, LSPRAY - a Lagrangian 
spray solver [14] and EUPDF - an Eulerian Monte 
Carlo PDF solver [15], were developed, thereby ex- 
tending our previous work on the Monte Carlo PDF 
and sprays [10] to unstructured grids as a part of 
NCC development. The combined unstructured 3D 
CFD/Spray/Monte-Carlo-PDF solver is designed to 
be massively parallel and accommodates the use of 
an unstructured mesh with mixed elements comprised 
of either triangular, quadrilateral, and/or tetrahedral 
type. 

A current status of the use of the parallel com- 
puting in turbulent reacting flows involving sprays, 
scalar Monte Carlo PDF, and unstructured grids was 
described in Ref. [11]. It outlined several numeri- 
cal techniques developed for overcoming some of the 
high computer CPU-time and memory-storage re- 
quirements associated with the use of Monte Carlo 
solution methods. The parallel performance of both 
the PDF and CFD modules was found to be excellent 
but the results were mixed for the spray computations 
showing reasonable performance on massively paral- 
lel computers like Cray T3D; but its performance was 
poor on the workstation clusters [10]. In order to im- 
prove the parallel performance of the spray module, 
two different domain decomposition strategies were 
developed and the results from both strategies were 
summarized [10-14]. 

Currently, most of the finite-difference tech- 
niques used for predicting the two-phase flows make 
use of the physics derived from the single- component 
liquid droplet studies with constant properties. How- 
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ever, it is well known that most of the gas-turbine fu- 
els are multicomponent mixtures of many compounds 
with a wide distillation curve [7,18]. The multicom- 
ponent nature of the liquid sprays is becoming evi- 
dent with the increasing need to use jet fuels derived 
from heavier petroleum compounds. The gasifica- 
tion behavior of a multicomponent fuel droplet may 
differ significantly over that of a pure single compo- 
nent fuel droplet [18]. Also, the calculation of the 
variable thermo-transport properties of the liquid- 
mixtures becomes more important at high pressures. 
The flame ignition characteristics such as the flame 
blow-off and extinction conditions could also be in- 
fluenced by the nonuniform concentration of the fu- 
els with different volatilities. However, the impor- 
tance of the multicomponent liquid fuels with vari- 
able properties received little attention in the mod- 
eling of comprehensive gas-turbine combustor spray 
computations. With this in mind, we have made the 
following improvements: 

(a) In order to deal with the gas turbine fu- 
els that are mixtures of many compounds, we have 
extended the spray formulation to multicomponent 
liquid sprays. This implementation also takes into 
account the effect of variable liquid properties. 

(b) In order to account for the nonideal gas be- 
havior under critical and supercritical conditions, we 
have completed the integration of the Peng- Robinson 
EOS into our CFD module and also, the implementa- 
tion of a high-pressure correction into the calculation 
of the transport properties in the gas phase. These 
modifications would enable the calculation of high- 
pressure flow properties in the gas-phase regardless 
of sprays. 

In this manual, we concentrate only on provid- 
ing the details of the spray module along with the 
high pressure corrections made to the gas-phase cal- 
culations. However, a discussion on the application 
of the joint scalar Monte Carlo PDF method to two- 
phase flows could be found elsewhere in [10-15]. Some 
of the salient features of the spray module are sum- 
marized below: 

• An efficient particle tracking algorithm was de- 
veloped and implemented into the Lagrangian 
spray solver in order to facilitate particle move- 
ment in an unstructured grid of mixed elements. 

• LSPRAY-II is currently coupled with an un- 
structured flow (CFD) solver of NCC [21-23], 
and an Eulerian-based Monte Carlo probabil- 
ity density function solver - EUPDF [15], which 


were selected to be as the integral components 
of the NCC cluster of modules. EUPDF was de- 
veloped for application with sprays, combustion, 
unstructured grids and parallel computing. 

• The spray solver receives the mean velocity and 
turbulence fields from the flow solver. The so- 
lution for the scalar (energy and species) fields 
could be provided by means of either a conven- 
tional CFD solver or a Monte Carlo PDF solver 
depending on the choice of the solver. 

• The spray solver supplies the spray source-term 
contributions arising from the exchanges of mass, 
momentum and energy with the liquid-phase 
to the flow solvers (CFD and/or Monte Carlo 
PDF). This information could be used in either 
conservative or non-conservative finite-difference 
formulations of the gas phase equations. 

The furnished code demonstrates the success- 
ful methods used for parallelization and coupling of 
the spray to the flow r code. First, complete details of 
the spray solution procedure is presented along with 
several other numerical issues related to the coupling 
between the CFD, LSPRAY-II, and EUPDF solvers. 
It is followed by a brief description of the combined 
parallel performance of the three solvers (CFD, EU- 
PDF, and LSPRAY-II) along with a brief summary 
of the validation cases. 

4 GOVERNING EQUATIONS FOR THE 
GAS PHASE 

Here, we summarize the conservation equations for 
the gas-phase in Eulerian coordinates [1]. These 
equations are valid for a dilute spray with a void frac- 
tion of the gas, 9 , close to unity. The void fraction is 
defined as the ratio of the equivalent volume of gas 
to a given volume of a gas and liquid mixture. This 
is done for the purpose of identifying the interphase 
source terms arising from the exchanges of mass, mo- 
mentum, and energy with the liquid-phase. 

The conservation of the mass leads to: 

[pVc],£ H“ [pVc^i] ,Xi — Smlc — ^ ^ Hk ^ k (1) 

k 

For mass conservation, the source term is given as 
a summation over different classes of droplets. Each 
class represents the average properties of a different 
poly disperse group of droplets. Here, represents 
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the number of droplets in a given class and rrik rep- 
resents the corresponding mass vaporization rate. 

For the conservation of the jth species, we have: 


\pVcyj\,t + [pVcDyj : Xi] ,Xi pV c Wj 


where 


$mls 


T, e i n k m k 

k 


( 2 ) 


^ Wj = 0 and ^ e j = 1 

3 3 

For the species conservation, the source term con- 
tains an additional variable, ej, which is defined as 
the fractional vaporization rate for species j. 

For the momentum conservation, we have: 


[pVc^i],t + \pV c UiUj\^ X j d" [pVcjjjCi 
[0VcTij] , X j — [(1 — $)Pc' 7 7ij] ,£j — Smlm = 

Z, n k m k u ki -^2 Pk r k n k u ki , t (3) 

k k 

where the shear stress rij in Eq. (3) is given by: 

7~ij = p[Ui : Xj + ^j,Xi\ 

For the momentum conservation, the first source 
term represents the momentum associated with liquid 
fuel vapor and the second represents the momentum 
change associated with droplet drag. 

For the energy conservation, we have: 


-[0V c p] t t = s m i e =^2 n k m k (h s - h,efi) (4) 

k 

Similarly, the energy conservation has a first source 
term associated with liquid fuel vapor and the second 
represents the heat loss associated with the latent 
heat of vaporization and it also contains an additional 
heat flux (loss or gain) to the droplet interior from the 
ambient. 

The main purpose of the spray solver is to cal- 
culate the source terms arising from the exchanges 
of the mass, momentum, and energy and, then, feed 


that appropriate information to the gas-phase solver. 
In the case of NCC, it suplies the source terms to the 
CFD solver and, also, to the Monte Carlo PDF solver 
if needed. 

The polynomial fits for the variable thermody- 
namic properties at low pressures are taken from the 
data set compiled by McBride et. el. [24]. The 
transport properties involving the thermal conduc- 
tivity and molecular viscosity for individual species 
is estimated based on the Chapman-Enskog colli- 
sion theory [25-27]. And Wilke’s formulae is used 
to determine the properties of mixture [25-27]. The 
binary-diffusion coefficients are determined based on 
the Chapman-Enskog theory and the Lennard-Jones 
potential [25-27]. 

5 HIGH PRESSURE EQUATION OF 
STATE 

In order to calculate the high pressure gas be- 
havior, the Peng-Robinson EOS (Equation- Of-St ate) 
is employed for a multi-component mixture in the fol- 
lowing form [16-17,28]: 


P = 


RT a m 

V^bZ~ U 2 + 25 mV -VL 


where 


(5) 


a m = E i E j ^yi(«i«i ) 1/2 ( 1 - hj), 
bm — Ej ■ 

l _ 0 . 07780 RT^ 
u i — p. 7 

* tc 

c = + / ‘“ <1 - T ‘P )]2 ’ 

Tir — T /Tic, 

f iuJ = 0.37464 + 1.54226^ - 0.26992^ 2 , 


tdi is known as the acentric factor of the molecules 
which is a measure of non-sphericity of the molecules, 
and kij is known as the binary interaction coefficient. 
However, the Peng-Robinson EOS is rewritten as a 
cubic EOS in terms of the compressibility factor, Z 
(= ^r), before it is solved: 

Z 3 -( 1 - B*)Z 2 + (A* - 2 R* - W* 2 )Z 

-A*B* + R* 2 + R* 3 = 0 (6) 
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where 


A* = and B* = b -f£. 

We have chosen the Peng-Robinson EOS be- 
cause of its simplicity and, more importantly, it 
proved to be very useful in the supercritical droplet 
vaporization studies of [19-20]. 

Table 1 lists various physical constants for some 
of the species that we found to be of interest in our 
spray computations. It contains values for the boiling 
temperature at normal pressure, the critical temper- 
ature, the critical pressure, the critical density, the 
latent heat of vaporization at normal pressure, the 
critical molar volume, the molar volume at normal 
pressure, and the acentric factor of the molecules. 
Most of this data is collected from [16-18] and is use- 
ful in both the evaluation of the Peng-Robinson EOS 
and the calculation of various other variable proper- 
ties. And Table 2 lists the binary parameters, kij, 
for a mixture of n-heptane, O 2 , IV 2 , CO 2 , & H 2 0 
used in a calculation involving the Peng-Robinson 
EOS. While the data for most of the binary pairs 
was obtained from various reference books, some of 
the missing data, however, is replaced with known 
data found for other binary pairs of molecules with 
similar molecular weights. 


and Q v = 1.0 for non-polar molecules. 

The effect of pressure on the thermal conduc- 
tivity of pure gases is determined by making use of 
the Stiel & Thodos method [16,19,30]: 


(A - \ n )YZl = 1.2210 _2 [ea;p(0.535p r ) - 1] (8) 

when p r < 0.5 

(A - \„)TZ b = 1.1410“ 2 [exp(0.67p r ) - 1.069] (9) 

when 0.5 < p r < 2.0 

(A - A„)r Z\ = 2.6010- 2 [ezp(1.155 / o r ) + 2.016] (10) 

when 2.0 < p r < 2.8 

where A is in W/(m.K), T = 210[ T ^ /3 ] 1/ ^ 6 , Z c is the 
critical compressibility, and p r is the reduced densitv 
p/Pc = V e /V. 

The effect of pressure & temperature on diffu- 
sion coefficients in gases is determined by making use 
of Takahashi correlation [16,19,31]: 


6 HIGH PRESSURE CORRECTIONS FOR 

GAS-PHASE TRANSPORT PROPERTIES . . = f(T r ,P r ) (11) 

{D A bP) + 


The effect of pressure on the viscosity of pure 
gases is determined by making use of the Reichenberg 
method [16,19,29]: 


— 1 + Q v 

P"n 


,3/2 


Ay Pi 

B v P r + (1 + Cy P, 




(7) 


where 


A v - 

= f-exp{a 2 T?) 

By 

— A v (/3i 

T r 

fa) 

C v - 

= %exp( 7 iT r c ) 

Dy 

= 7j±exp(\?.Td) 

ai = 

= 1.982510 -03 

ol 2 — 

5.2683 

a 

= -0.5767 

A = 

= 1.6552 

P2 = 

1.2760 



7i = 

= 0.1319 

72 = 

3.7035 

c 

= -79.8678 

Ai = 

= 2.9496 

a 2 = 

2.9190 

d 

= -16.6169 


where Dab — diffusion coefficient, cm 2 /s, P = 
pressure, bar, the superscript + indicates that low- 
pressure values are to be used, and the reduced pres- 
sure and temperature in the above equation are cal- 
culated as follows: 


T — X 

± r — j' 

T c = VaT c a + VbT c b 

P — — 

±r p c 

Pc = VaPcA + VbPcB 


This correlation is valid for a system involving 
a trace solute diffusing in a supercritical fluid. It is 
shown to provide satisfactory results but it is based 
on a very limited set of available experimental data. 
The graphical representation of the Takahashi func- 
tion is given in Fig. 1. 
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Table 1. Physical constants. 


Species 

T n b 

(K) 

T c 

(K) 

p c 

(atm) 

Pc 

(Kg/m 3 ) 

Ikin 

(KJ/kg) 

V e ^ 

(cm 3 /g-mole) 

V n 

(cm 3 /g-mole) 

0J 

a 

(A °) 

a / n 

(K) 

C 0 Hu 

341.9 

507.4 

30.0 

660.0 

334.8 

370.0 

140.06 

0.296 

5.949 

399.3 

c 7 h u 

371.6 

540.2 

27.0 

682.0 

316.3 

432.0 

162.00 

0.351 

6.297 

419.031 

C$His 

398.82 

568.8 

24.6 

718.5 

301.3 

492.0 

188.8 

0.394 

6.62 

488.15 

ClO-^22 

447.3 

617.6 

20.8 

728.3 

276.1 

603.0 

233.68 

0.490 

7.16 

540.06 

Cl 2 #26 

489.5 

658.3 

18.0 

748.0 

256.3 

713.0 

278.54 

0.562 

7.655 

583.68 

C U H W 

526.7 

694.0 

16.0 

763.0 

240.1 

830.0 

326.62 

0.679 

8.067 

629.08 

n 2 

77.4 

126.2 

33.9 

807.1 

197.6 

90.1 

31.87 

0.039 

3.681 

91.5 

o 2 

90.2 

154.6 

50.4 

1135.7 

212.3 

74.4 

26.08 

0.025 

3.433 

113.0 

co 2 

00.0 

304.1 

73.8 

000.0 

000.0 

94.0 

33.32 

0.239 

3.996 

190.0 

h 2 o 

373.2 

647.3 

221.2 

958.1 

2257.2 

57.1 

19.76 

0.344 

2.641 

809.1 


Table 2. Binary parameter constants, kij. 


C 7 H i 6 
0=1) 

o 2 

(j=2) 

n 2 

0=3) 

C0 2 

0=4) 

h 2 o 

0=5) 

c 7 h 16 

(i=l) 

0.0000 

0.1321 

0.1440 

0.1000 

0.1484 

o 2 

(i=2) 

0.1321 

0.0000 

-0.0119 

-0.0289 

0.0910 

N . 2 
(i— 3) 

0.1440 

-0.0119 

0.0000 

-0.0170 

0.1030 

CO 2 
(i=4) 

0.1000 

-0.0289 

-0.0170 

0.0000 

0.1200 

h 2 o 

(i=5) 

0.1484 

0.0910 

0.1030 

0.1200 

0.0000 
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P/(D ab P) 


as given by Eq. (46). The droplet Reynolds number 
is given by: 



0.0 1.0 2.0 3.0 4.0 5.0 6.0 


Reduced pressure, P r 

Figure 1 Takahashi correlation showing the varia- 
tion of the binary diffusion coefficient versus reduced 
pressure at different reduced temperatures. 


Re k = 2 r -^- [(u 3 +u' g - u k ) . (u g +u g - u*)] 

P-gs 

(15) 

By following the approach taken from KIVA 
[33], a gas turbulence velocity, v! g , is added to the 
local mean gas velocity when calculating the droplet 
drag and vaporization rates. The gas velocity fluctu- 
ations is calculated by randomly sampling a Gaussian 
distribution with mean square deviation, 2/3 k. The 
Gaussian is given by 


G{u' g ) = (4/37rfe) 3 / 2 exp[-3\u' g \ 2 /Ak\ (16) 

The gas fluctuating component is calculated once at 
every turbulence interaction time, ttur, and is oth- 
erwise held constant [33]. The correlation time is 
taken to be the minimum of either the eddy time or 
the transit time taken by the droplet to traverse the 
eddy. It is given by 


7 LIQUID-PHASE EQUATIONS 

Here, we summarize the governing equations for 
the liquid-phase based on a Lagrangian formulation 
where the equations for particle position and velocity 
are described by a set of ordinary differential equa- 
tions. For the particle position of the kth droplet 
group, we have: 


dxj k 

dt 


= U-ik 


For the droplet velocity: 


( 12 ) 


t tur — 7TIXTI I , Cft, 


k fc 3 / 2 

e e \u g + u' g -Uk\- 


(!• 


— ) 

uL — Uh 1/ 


(17) 


where cu is an empirical constant with a value of 


0.16432. 

The liquid mixture density, pk, in Eq. (13) is 
calculated as follows [16,18]: 


pk — ^^VkiPki ( 18 ) 

i 

and the individual component liquid density is given 
by: 


du-ik 

dt 


3 CpptgsRek 
16 purl 


\u ig + u' g 


Wife] 


where 


( 13 ) _ M ki 

Pki — T/ 

F/ci 

where the molar volume, T4i, is 




(14) 

and 


V ki = Va (0.29056 - 0.8775o/.;) c '’ 


(19) 


( 20 ) 


According to Yuen and Chen [32] , the droplet drag for g 

a vaporizing droplet could be calculated by a solid- c v = [l — yp- 

sphere drag correlation but suggested using a correc- Cl 

tion in the evaluation of the droplet Reynolds number The droplet regression rate is determined from 

for the average viscosity based on a 1/3 weighting rule one of three different correlations depending upon the 
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droplet Reynolds number range. The first correlation 
as given by Eq. (21) is based on a gas-film analysis 
developed by Tong & Sirignano [34] . It is based on a 
a combination of stagnation and flat-plate boundary- 
layer analysis and is valid for Reynolds numbers in 
the intermediate range. The last two correlations as 
given by Eqs. (22) & (23) are taken from Clift et 
al [35] and are valid when Re k < 20. In fact, when 
the droplet Reynolds number goes to zero, Eq. (23) 
becomes identical to the droplet regression rate of a 
vaporizing droplet in a quiescent medium [27]. 


ds k 

dt 



r 2 l 1 / 2 
-Re k f(B k ) 

7 r 


if Re k > 20 


( 21 ) 


ds k 

dt 


= [l + (1 + Re k ) 1/3 ] Re° k 077 ln( 1 + B k ) 

^ ’ (22) 


if 1 < Re k < 20 


ds k 

dt 


PgDg 

Pk 


1 + (1 + iJe*;) 1 / 3 ] ln(l + B k ) 


(23) 


if Re k < 1 


where B k is the Spalding transfer number and is given 
by: 


B k 


(Vfs ~ Vf) 
(1 - Vfs) 


(24) 


and yf s is given as a summation over the fuel-species 
mass fractions at the droplet interface: 


yfs = y ™ ( 25 ) 

i 

and yf is given as a summation over the fuel-species 
mass fractions in the ambient: 


w = &/i ( 26 ) 

i 

and the function f(B k ) is similar to that of a Blassius 
function [1, 36] and is obtained from an analysis simi- 
lar to that of Emmon’s boundary-layer flow over a flat 
plate with blowing [36]. The range of validity of this 
function was extended in Raju and Sirignano [7] to 


consider the effects associated with a boundary-layer 
flow with suction. 

The internal droplet temperature is determined 
based on a vortex model [34]. The governing equation 
for the internal droplet temperature is given by: 


9T k _ A ! 
dt C v ipir\ 


a- 


d 2 T k 
da 2 


+ (1 + C(t)a) 


m 

da 


(27) 


where 


C{t) = 


3_ 

17 


Cpipi 

A/ 



(28) 


where a represents the coordinate normal to the 
stream-surface of a Hill’s Vortex in the circulating 
fluid, and C(t) represents a nondimensional form of 
the droplet regression rate. The initial and boundary 
conditions for Eq. (27) are given by: 


t — t 


injection •> 


T k = T kj0 


a = 0, 


m 

da 


17 


Cpipi 1 2 dT k 

A , J T ' k dt 


(29) 

(30) 


a — 1, 


m 

da 


3 Pk r J J 1 ds k 

~32T l [lk ' €ff ~ lk] lT 


(31) 


where a = 0 refers to the vortex center, and a = 1 
refers to the droplet surface, and the mixture latent 
heat of vaporization, l k , is given by 


h — ^2 Ci ^ ki (32) 

i 

and the individual component latent heat of vapor- 
ization, l k i , is given by [16,18]: 


hi = hin(j^ 


-?fcX0-38 

-tJ 


(33) 


and the droplet boiling temperature is given by 


rp _ l kin Mi / R U /q 

H ~ IkinMf/iRuh^-lniP) K } 

and, finally, the effective latent heat of vaporization, 
l ki eff, is defined as: 

h,eff=h+ 47T^(^) (35) 

m k \ dr J s 

It is a very useful parameter as it represents the total 
energy loss associated with the latent heat of vapor- 
ization in addition to the the heat loss to the droplet 
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interior. h,eff is calculated by means of the following 
relationship [18]: 


/ _ Cp{T g Tks) 

k ' eff ~ (1 + Bk) l t Le — 1 


(36) 


Similar to the internal thermal transport, the 
internal mass transport of a multi- component fuel is 
given by: 


&yki ^7 Dk 

~dT ~ ~rf 


d 2 Vki , , n(±\~\ 

ol- r-r- + (1 + C (t)a) 


da 2 


da 


(37; 


The initial and boundary conditions for Eq. (37) are 
given by: 


^ t injection-) Vki Vki,o 


0 = 0 , 


dyki 

~d^~ ~ 17 


r -2 -1 


D k . 


dyki 

dt 


(38) 

(39) 




— Vis + (l — y/s ) 


Vis - yfi 
yfs - yf 


(45) 


It is noteworthy that the thermodynamic and 
transport properties at the gas film are calculated at 
the temperature and composition as determined by 
the following one- third rule: 


1 2 

< fiavg — 3 ^ (46) 

The correlations for the gas-phase thermody- 
namic and transport properties are described in Sec- 
tions 4 & 6 . In a similar way, the liquid-phase ther- 
modynamic and transport properties are determined 
based on the correlations described in the next sec- 
tion. 


8 LOW PRESSURE LIQUID 
THERMODYNAMIC & TRANSPORT 
PROPERTIES 


a = 1 , 



3 1 r 1 dSk 

~32 Tr [Vkis - Cii ~dt 


The specific heat at constant pressure, C p i, 
(40) thermal conductivity, A/, and viscosity, pi, are evalu- 
ated by means of the following expressions: 


By knowing the mass fractions of the liquid 
species at the droplet surface, the corresponding mole 
fractions are determined by 


yiks/Mi 

XikS ~ Ei Vyks/Mi 


(41) 


At the droplet interface, the mole fractions of 
the gas species are obtained by means of Raoult’s 
law: 


%is — p%iksPis (42) 

where the partial pressure, Pi S , is determined by 
means of the Clausius-Clapeyron relationship: 


P is = exp 


Iki j 


V 

R u 

\T bi 

Tks) . 


(43) 


Then, the corresponding mass fractions in the gas- 
phase at the droplet interface are given by 


yis — 


XisMj 

Af a (l — %is) 4” E?, MiXis 


(44) 


where M a is the molecular weight of the gas excluding 
fuel vapor. And the fractional mass vaporization rate 
of liquid species, e$, is given by 


Cpi — Cpio + C p i\T -f- CpioT 2 + Cpt^T 3 + C p iaT a 

(47) 


Xi = Xio + X n T + A i2 T 2 + A/ 3 T 3 + A; 4 T 4 + A / 5 T 5 

(48) 


In in — nio P P11/T + pizT + P13T 2 (49) 


where C p i is in J/(kg K), pi in (pPk s), and A / in 
(W/m K). 

Tables 3-5 provide the polynomial constants 
used in Eqs. (47) -(49) for some of the species that 
we found to be of interest in our spray computations. 
Table 3 provides the constants for the liquid specific 
heat, Cpi , Table 4 for the liquid thermal conductivity, 
A/, and Table 5 for the liquid molecular viscosity, pi. 
These tables are compiled with the data taken mostly 
from the references of [16,18]. 

The binary diffusion coefficient, Dij, is evalu- 
ated as follows [16]: 


Da = 


K dif T 

NVP 


(50) 
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Table 3. Polynomial constants for liquid specific heat. 

Fuel 

CplO 

Cpll 

Cpl2 

CplZ 

CplA 

c 6 h 14 

2.4169 

-5.9866e-03 

2.0959e-05 

-8.4489e-09 

0.0 

c 7 h 16 

4.8227 

-3.6980e-02 

1.6777e-04 

-3.0987e-07 

2.2081e-10 

C»H is 

9.2189 

-8.8314e-02 

3.8869e-04 

-7.2539e-07 

5.0776e-10 

Cl()H22 

4.7991 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 

^12^26 

4.7900 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 

Cl 4 i/30 

4.7991 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 


Table 4. Polynomial constants for liquid thermal conductivity. 

Fuel 

Ajo 

Aii 

A/2 

A/3 

A/4 

A/5 

c 6 h 14 

0.37078 

-5.4313e-03 

4.628e-05 

-1.8002e-07 

3.2243e-10 

-2.1832e-13 

C r H lft 

0.13236 

9.4441e-04 

-6.588d-06 

1.4617e-08 

-1.1244e-ll 

0.0 

C S H i 8 

0.25652 

-7.5401e-04 

1.5872e-06 

-1.6795e-09 

-1.3375e-16 

0.0 

C 10 H 22 

0.22179 

-2.3699e-04 

-6.94e-07 

2.0415e-09 

-1.5741e-12 

0.0 

C 12 H 26 

0.17609 

4.2463e-05 

-7.4467e-07 

6.9446e-10 

0.0 

0.0 

CuH 3 0 

0.18801 

-9.1399e-05 

-2.1464e-07 

1.1655e-10 

0.0 

0.0 

n 2 

-2.629E-1 

-1.545E-3 

-9.450E-7 

0.0 

0.0 

0.0 

O 2 

2.444E-1 

-8.813E-4 

-2.023E-6 

0.0 

0.0 

0.0 

CO 2 

4.070E-1 

-8.438E-4 

-9.626E-7 

0.0 

0.0 

0.0 

h 2 o 

-3.838E-1 

5.254E-3 

-6.369E-6 

0.0 

0.0 

0.0 


Table 5. 

Polynomial constants for liquid molecular viscosity. 

Fuel 

Mzo 

m 

^12 

M/3 

C 6 H 1A 

-4.034E+00 

8.354E+02 

0.0000000 

0.0000000 

c 7 h u 

-4.325E+00 

1.006E+03 

0.0000000 

0.0000000 

C S H i 8 

-4.333E+00 

1.091E+03 

0.0000000 

0.0000000 

CioH‘22 

-4.460E+00 

1.286E+03 

0.0000000 

0.0000000 

C 12 H 26 

-4.562E+00 

1.454E+03 

0.0000000 

0.0000000 

C 14 H 5 0 

-4.615E+00 

1.588E+03 

0.0000000 

0.0000000 

n 2 

-2.795E+01 

8.660E+02 

2.763E-01 

-1.084E-03 

0 2 

-4.771E+00 

2.146E+02 

1.389E+02 

-6.255E-05 

CO 2 

-3.097E+00 

4.886E+01 

2.381E-02 

-7.840E-05 

h 2 o 

-2.471E+01 

4.209E+03 

4.527E-02 

-3.376E-05 
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where K^if = 8.210 _8 (1 + [^] 2 / 3 )- One should be 
careful in using this approximation as it is based on 
a scarce set of available experimental data. 

The specific heat for a multicomponent mixture 
is given by: 

n 

Cpm ~ ^ ^ Ui^pi ( 51 ) 

i=l 

And the thermal conductivity of a multicomponent 
mixture is calculated by means of the Li method [16]. 

n n 

^ m — ^ ^ ^ ^ ( 52 ) 

i— 1 j = 1 

where 

Xij = 2(A- 1 + A- 1 ) _1 

A. - 


where Xi is the mole fraction of the species i, fa is a 
volume fraction of the ith species, and Vi is the molar 
volume of the pure fluid. 

9 DETAILS OF DROPLET FUEL 
INJECTION 

The success of any spray model depends a great 
deal on the specification of the appropriate injector 
exit conditions. However, a discussion involving the 
physics of liquid atomization is beyond the scope of 
this subject matter. A great deal of research still 
needs to be done before we would be able to model 
the underlying atomization characteristics of differ- 
ent injectors from the first principles. However, in 
our present computations, we rely on some widely- 
used correlations to provide for the initial spray con- 
ditions. And the liquid fuel injection is simulated by 
introducing a discretized parcel of liquid mass in the 
form of spherical droplets at the beginning of every 
fuel-injection time step, A tu. 

The spray computations facilitate the use of 
multiple fuel injectors. The same or a different type 
of liquid fuel can be specified for each one of different 
injectors. The initial droplet temperature is assumed 
to be the same for all different droplet groups of any 
given injector. The following three variables play an 
important role in simulating the injector initial con- 
ditions: 

• no_of-holes() - the number of holes per injector, 


• no_of_streams() - the number of droplet streams 
per hole - it is introduced to distinguish the 
initial velocity variation within different droplet 
classes arising from the geometric considerations 
of a chosen spray, & 

• no_of_droplet .groups () - the number of droplet 
groups per stream - it is introduced to distin- 
guish the droplet-size variation within different 
droplet classes of a polydisperse spray. 

The total number of droplet groups introduced at the 
beginning of every different injection time step for 
a given injector is therefore equal to a value given 
by the multiplication of these three variables. Fur- 
ther details on some of the spray input variables 
can be found in the spray input file, nccJnjector.in.l 
of Table 6, where the integer number followed by 
the last dot represents the injector number, which 
in this case happens to be one. The table pro- 
vides a complete description of the following in- 
put variables: out_string, tdrop(n_i), ymki(n_i,n_l), 
spray -table, steady_spray .model, no_of_holes(nd), 
no_of_streams(n_i), no_of_droplet_groups(n_i), lmdis, 
smdm, cone, ((xdnj(nx), y_inj(nx), z_inj(nx), 
flowf(nx), v_inj(nx), alphadnj (nx) , beta_inj(nx), 
thetadnj(nx), dtheta_inj(nx), swlr_angle(nx), 
size_min(nx), size _max (nx)), nx=l,no_of_holes(n_i)). 

The initial droplet distribution for each injec- 
tor could be defined either by a complete specifica- 
tion of the initial conditions through a spray table or 
through the specification of some input parameters 
that invoke certain pre-defined correlations. 

When a spray table is defined, one should pro- 
vide the information on the (x,y z) components of 
the initial droplet location, the (u, v, w) compo- 
nents of the droplet velocity, and the initial mass flow 
rate associated with each different droplet group as 
described in ncc_spray_table.l of Table 7. This ta- 
ble provides a complete description of the following 
variables: nos(n_i), (ni,xx_inj, yy_inj, zzdnj, uu_inj, 
vv_inj, ww_inj, r_inj, fld_d), ,ni=l,nos(n_i)). The ini- 
tial inputs specified through a spray table should be 
representative of the integrated averages of the ex- 
perimental conditions [10-13]. 

The remainder of the discussion in this section 
applies to the case when a spray table is not de- 
fined. In which case, we need to specify several pa- 
rameters including no_of_holes(), no_of_streams(), & 
no_of_droplet_groups(). And depending on what is 
specified for the input of the logical variable, lmdis, in 
ncc_injector.in.l, the droplet-size distribution within 
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Table 6. Typical nccJiquid_injector.in.01 file. 

Input file content 

comments 

heading 

title of a description of this file 

heading 

title of property 

nolc(nJ) 

denotes the total number of initial liquid components 
in the nJ-th injector. 

heading 

title of property 

out_string 

chemical symbols of initial liquid species: e.g., C6H12. 

(ymki(nJ,n_l), 

nJ=l,nolc(nJ)) 

initial mass fractions of liquid species 

heading 

title of controlling parameters 

tdrop(nJ) 

initial droplet temperature 

heading 

title of controlling parameters 

spray_table, 
steady_spray .model 

If spray-table = .true., initial droplet location, velocity, size, 
and mass flow rate are defined in the file - nccJiquid-table.in.l. 
Otherwise they are determined based on some form of specified 
spray correlations and configurations. 

If steady .spray .model = .true., it invokes a steady state 
spray model commonly used in many spray codes, where 
whenever the spray solver is called, after first introducing 
a new group of spray particles, it continues with the 
liquid-phase computation until after all the particles are taken 
out of the computational domain. (NOTE: It is NOT 
recommended to use this option as a steady state calculation 
could be better arrived at by making use of some features of 
the unsteady option.) The solution from the unsteady option 
(steady .spray .model = .false.) is determined based on the 
values assigned to the controlling time steps - dtml, dtil, 

& dtgl, which are internal to the spray solver. 

heading 

title of controlling parameters 

no-of-holes(nJ), 
no_of_streams(nJ), & 
no_of_droplet -groups (ni) 

no_of_holes(nJ) = The number of holes per injector. 

no_of_streams(nJ) = The number of droplet streams per hole. 
This variable is introduced mainly to distinguish the variation 
in the droplet groups based on angular orientation. 

no_of_droplet-groups(nJ) = The number of droplet groups per 
stream. This variable is introduced mainly to distinguish the 
droplet groups based on the droplet-size variation. 

NOTE: when spray -table = .true., both no_of_holes and 
no_of_streams are set equal to 1 and all the different droplet 
groups are grouped into a single variable, no_of_droplet_groups(). 
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Table 6. Typical nccJiquid_injector.in.01 file (continued). 

Input file content 

comments 

heading 

title of controlling parameters 

lmdis, smdm, & cone 

If lmdis = .true., it invokes a correlation for droplet size distribution. 
Otherwise it assumes a uniform droplet distribution between the maximum 
and minimum initial droplet sizes as defined by size_max() & size_min(). 

smdm = Sauter mean diameter 

If cone = .true., it activates a 3D solid or hollow cone spray configuration 
as shown in Fig. 2. Otherwise it activates a 2D configuration of either an 
axis-of-symmetry case of Fig. 3 or a planar case of Fig. 4 depending on the 
value of the logical variable - axisymmetric. 

heading 

title of controlling parameters 

( (x-inj (nx) ,y _inj (nx) , 
z Jnj (nx) ,flowf(nx) , 
v Jnj (nx) , alpha Jnj (nx) , 
beta Jnj (nx) , theta Jnj (nx) , 
dthetaJnj(nx), 
s wlr .angle (nx), 
size_min(nx),size_max(nx)), 
nx= 1 ,no_of Jioles (n J) ) 

(xJnj(nx),yJnj(nx),zJnj(nx)) = spatial coordinates of the initial injector 
location. 

flowf(nx) — injector mass flow rate, (units - kgm/s for 3D & axis-of- 
symmetry & kgm/s/m for 2D planar. REMEMBER FOR AXIS-OF- 
SYMMETRY, IT IS THE TOTAL FLOW RATE OVER 360 DEGS.) 

alphaJnj(nx) = angle of rotation towards z-axis. 

betaJnj(nx) = angle of rotation towards y-axis. 

theta Jnj (nx) = cone angle. 

dtheta Jnj (nx) = half-cone angle. (NOTE: Although dtheta Jnj (nx) 

= thetaJnj(nx)/2 for SOLID CONE SPRAY, it is invoked by 
setting dtheta Jnj (nx) either equal to 0 or thetaJnj(nx)/2.) Refer 
to Figs. 2-4, for a better understanding of the angular representation. 

swlr .angle (nx) = swirl angle. 

size_min(nx) & size_max(nx) = The variables are associated with 
the droplet size distribution of lmdis = .false. 
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Fig 2a. Geometrical details of fuel injection for 
a 3D solid or hollow cone spray. 



(j) = azimutal angle, deg 
number of rays = 8 (<|)=45) 

For a 3D cone, 
no_of_str earns ( ) should 
defined as a multiple of 8. 
number along each ray = 
no_of_st reams 0/8. 


Fig 2b. Initial spray particle orientation in a 
circular corss-section . 
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0 = cone angle, deg 
d0 = hal f-cone angle 



Fig 3a. Geometrical details of fuel injection for 
an axis-of -symmetry case. 


f 00 

I 

f For axis of symmetry: 

f no_of_streams ( ) is distributed 

• 01 between 01 to 00 . 


I 

: C 


Fig 3b. Initial spray particle concentration in an 
axis-of -symmetry case. 
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Table 7. ncc_spray_table.in.01 file. 

Input file content 

comments 

nos(nd) 

denotes the total number of droplet groups in the n_i-th 
injector. 

(ni,xx_inj,yyJnj, 
zz_inj ,uu_inj ,vvdnj , 
wwdnj ,r_inj ,fld_d) 
,ni=l,nos(n_i)) 

ni = number of the droplet group. Its value ranges between 
1 to nos(n_i). 

(xxdnj,yydnj,zzdnj) = spatial coordinates. 

(uu_inj,vv_inj, wwdnj) = velocity components, 
rdnj = droplet size in radius. 

fld_d = mass flow rate of the ni-th droplet group. (NOTE: 
SUMMATION OF fld.d OVER ALL THE nos(nd) 
DROPLET GROUPS IS EQUAL TO THE TOTAL MASS 
FLOW RATE OF THE INJECTOR.) 


each one of the streams is determined based on either 
one of the two options: 

• In one option, it is calculated based on a corre- 
lation typical of those widely used in describing 
the initial droplet size distribution [4]: 


— = 4.21 x 10 6 
n 


^32 


3.5 


— 16.98 ( 3 ” ) ' dd 
dyi 

(53) 


where n is the total number of droplets and dn is 
the number of droplets in the size range between 
d and d + dd. This correlation also requires the 
specification of Sauter mean diameter, c? 32 . Fig. 
5 shows the droplet size distribution generated 
by Eq. (53) for a case studied in [10]. The solid 
line shows the droplet number variation versus 
drop size and the dashed line shows the inte- 
grated mass variation with drop size. The drop 
size distribution within the spray is represented 
by a finite number of droplet classes as given by 
the variable, no_of .droplet -groups (). 



Figure 5 Droplet-size distribution. 


• In the second option, the initial droplet sizes are 
distributed evenly between two specified maxi- 
mum and minimum drop sizes - size_max() & 
size_min() as specified in Table 6. And the size 
interval depends on the value specified for the 
variable, no_of_droplet_groups(). However, this 
option is applicable only in certain special cases. 


Depending on what is specified for the logical 
variable, cone, of Table 6, the droplet velocity dis- 
tribution amongst various streams of a given hole is 
calculated by assuming the spray to be either a solid 
or hollow cone spray. A graphical illustration of three 
different cone configurations are shown in Figs. 2 to 
4. Figs. 2a & 2b refer to a 3D case, Figs. 3a &; 3b 
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to an axisymmetric case, and Figs. 4a & 4b to a 2D 
planar case. It is noteworthy that in an axisymmet- 
ric case, the x-axis is assumed to be aligned with the 
axis-of-symmetry. 

Fig. 2a shows the geometric details of a hollow 
cone spray in 3D where 6 is known as the cone angle, 
d6 is the half-cone angle, and for a solid cone spray 
dO = 0/2. And a represents the angular rotation to- 
wards the z-axis and /? is the angular rotation towards 
the y-axis. Fig. 2b shows initial spray stream orien- 
tation in a circular cross section. We try to simulate 
the spray by a finite number of streams as given by 
the variable, no_of_streams(). Each one of the dark 
circles in Fig. 2b represents a different stream. The 
streams are distributed evenly along each one of the 
different rays which are separated from each other 
with an angle of separation as given by <f>. Currently, 
we have hard-coded the angle of separation to be 45° 
which means we have restricted the number of rays 
to be eight in 3D. Therefore, when specifying a num- 
ber for no_of_streams(), it should be borne in mind 
that it should be a multiple of eight. And the num- 
ber of streams along each one of the rays, therefore, 
becomes equal to no_oLstreams()/8. In Fig. 2b, the 
no_of_streams() has a value of 32 and, therefore, the 
number of streams along each one of the rays for this 
case is 4. 

Fig. 3a. shows the geometric details for an 
axisymmetric case. Since the computations are per- 
formed only in the first quadrant of the x-y plane, all 
of the specified number of streams, no.oLstreamsQ, 
are distributed evenly between the lower and upper 
halfs as shown in Fig. 3a and 3b. Here, the upper half 
refers to the region between 01 to 00 and the lower 
half refers to LO to LI. It is noteworthy that each 
droplet group in a given stream represents a circular 
ring of liquid for an axisymmetric case. 

The geometric details for a 2D case are shown 
in Figs. 4a & 4b. Here, it is noteworthy that 
each droplet group in a given stream represents a 
planar sheet of liquid. All the specified number of 
streams, no_of_streams(), are distributed evenly over 
both sides of the cone center as shown in Fig. 4b. 

For each one of the different injector holes, 
no_of_holes(), it requires the specification of the fol- 
lowing parameters as described in ncc-injector.in.l 
of Table 6 (However, because of the geometric dif- 
ferences that exist between 3D, 2D, and axisymmet- 
ric configurations, some of the input parameters may 
have different units as noted below): 


• The initial (x, y, z) coordinates of the hole loca- 
tion. 

• The mass flow rate per hole - however, the defi- 
nition of the units for the injector mass flow rate 
per hole differs: it is kgm/s for 3D & the axis-of- 
symmetry and kgm/s/m for 2D planar (NOTE: 
In an axis-of-symmetry, the specified mass flow 
rate per hole refers to the entire mass flow rate 
over 360 degrees). 

• The following variables define the angular orien- 
tation: Oii n j = angle of rotation towards z-axis, 
fiinj = angle of rotation towards y-axis, 0 in j = 
cone angle, & dOi n j = half-cone angle (NOTE: 
Although dO in j = 0 i n j / 2 for a solid cone spray, 
a specified value of zero for dd{ n j also invokes a 
solid-cone spray configuration). 

• The variable, swlr_angle(), allows a means to 
specify the tangential component in the case of 
both 3D and axisymmetric sprays. 

10 SPRAY SOLUTION ALGORITHM 

In order to evaluate the initial conditions that 
are needed in the integration of the liquid-phase equa- 
tions, we first need to know the surrounding gas- 
phase properties at each particle location. But in or- 
der to evaluate the gas-phase properties, it is first nec- 
essary to identify the computational cell in which a 
given particle is located. It is a trivial task to track a 
particle in the regular rectangular coordinates. How- 
ever, the particle tracking becomes complicated when 
the computational cell is no longer rectangular in the 
physical domain and it becomes even more compli- 
cated when the particle search is performed within 
the context of parallel computing. 

We have developed and implemented an effi- 
cient particle-tracking algorithm for use with parallel 
computing in an unstructured grid. The search is ini- 
tiated in the form of a local search originating from 
the computational cell in which the same particle was 
found to be located in the previous time step. The 
location of the computational cell is then determined 
by first evaluating the dot product of x pc . a n = \x pc \ 
|a n | cos ), where x pc is the vector defined by the 
distance between the particle location and the center 
of the n-face of the computational cell, a n is the out- 
ward area normal of the n-face as shown in Fig. 6, 
and <fi is the angle between the two vectors. 
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A simple test for the particle location requires 
that the dot product be negative over each and ev- 
ery one of the n-faces of the computational cell. If 
the test fails, the particle search is carried over to the 
adjacent cells of those faces for which the dot prod- 
uct turns out to be positive. Some of those n-faces 
might represent the boundaries of the computational 
domain while the others represent the interfaces be- 
tween two adjoining interior cells. The search is first 
carried over to the adjacent interior cells in the di- 
rection pointed out by the positive sign of the dot 
products. The boundary conditions are only imple- 
mented after making sure that all other remaining 
possibilities point towards a search exterior of the 
computational domain. This implementation ensures 
against any inadvertent application of the boundary 
conditions before the correct interior cell could be 
identified. 



Fig. 6 A vector illustration used in the particle 
search analysis. 


After the gas-phase properties at the particle 
location are known, the solution for the ordinary dif- 
ferential equations of particle position, size, and ve- 
locity are advanced by making use of a second-order 
accurate Runge-Kutta method. The partial differen- 
tial equations governing the droplet internal thermal 
and mass transport are integrated by making use of 
a fully implicit Newton-Raphson iteration method. 

Finally, the liquid-phase source terms of the gas- 
phase conservation equations (1-4) are evaluated by 
making use of a time-averaging method. 


11 THE FLOW STRUCTURE OF THE 
SPRAY CODE & A DESCRIPTION OF 
THE TIME-AVERAGING SCHEME USED 
IN THE CALCULATION OF THE 
SOURCE TERMS 

• In order to know more about the time-averaging 
method, we need to know first about the three 
different time steps that are internal to the spray 
code: A t mh At«, and A t gl . 

A t ml - the actual time step used in integrating 
the liquid-phase equations which is determined 
based on the smallest of the different time scales 
associated with various rate-controlling phenom- 
ena of a rapidly vaporizing droplet. Some of 
the limiting scales axe the average droplet life- 
time, the time it takes for the droplet to traverse 
the local grid spacing, the time it takes for the 
droplet internal temperature to reach the liquid 
boiling temperature, & a relaxation time scale 
associated with droplet drag. The time-scale re- 
striction based on these criteria can become quite 
severe for smaller drops - for drops of sizes less 
than a micron. 

A tu ~ the injection time step. It is the time 
step at which a new discretized parcel of different 
droplet groups are introduced into the computa- 
tion. 

A tgi - the global time step. Its introduction 
seems to provide better convergence in both un- 
steady and steady-state computations. 

• When the spray solver is called it advances the 
liquid phase equations over a number of itera- 
tions as determined by the ratio of A t g i/ At m i . 

• It then evaluates the time-averaged contribution 
of the liquid-phase source terms, S g i , of the gas- 
phase governing equations (1-4) as follows: 

M 

S 3 i=E 

771=1 

where 


M 

^ A t m i = (55) 

m=l 

• The values for A t mh A t gh & A t u are specified 
in the input file, ncc Jiquid-solver.in, of Table 8. 


At ml 


At 


& 


Sml 


(54) 
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Fig. 7 The flow structure of LSPRAY-II. 
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• In steady-state computations, it is recommended 
to use for both A t g i and A tu a value of about 
1 ms which is roughly equivalent to the aver- 
age lifetime of the droplets for a typical reacting 
spray encountered in conventional low-pressure 
gas-turbine combustors. 

The averaging scheme could be explained better 
through the use of a flow chart shown in Fig. 7. 
The main spray solver is invoked with a controlling 
routine, DCLR, which, then, executes the following 
steps: 

L It first initializes the source terms to zero. 

2. Updates the global time, t g t , based on At^. 

3. Checks to see if t m i < t g t < t m i + A t m i. If it 
is, it returns control over to the calling routine 
and supplies the other flow solvers, e.g., flow or 
EUPDF, with the source terms, S g i , of Eqs. (1)- 
(4). If not, it proceeds with the next step. 

4. Checks to see if it it is time to introduce a new 
group of particles. 

5. Proceeds with solving the liquid-phase equations 
with calls to the following routines: 

• Calls the particle tracking routine and as- 
signs particles based on the parallel strategy 
implemented. 

• Interpolates gas-phase properties at the 
particle location. 

• Advances liquid-phase equations and, then, 
deletes any particles that are no longer 
needed in the computations. 

6. Evaluates the liquid-phase source-term contribu- 
tions, S m j, of Eq. 54. 

7. Updates the time, t mi f, based on A £ mf . 

8. It then goes back to step (3) and repeats 

the whole process again until the computa- 
tions are completed over a global time step: 
X) m =i = A t g i. 


12 IMPLEMENTATION OF THE 
BOUNDARY CONDITIONS IN LSPRAY-II 

The spray code supports most of the boundary 
conditions that are in use in the current version of 
the NCC CFD module but not all. Here, we would 
like to highlight some of the boundary conditions: 

• In case of the droplet impingement with the com- 
bustor walls, the droplets may evaporate, move 
along the wall surfaces, and/or reflect with re- 
duced momentum. The physics of droplet im- 
pingement with the walls is not completely un- 
derstood. In our present calculations, it is as- 
sumed that the droplets, after having lost most 
of their momentum upon impingement with the 
walls, move along the wall surfaces with a veloc- 
ity equal to that of the surrounding gas. 

• The implementation of the periodic boundary 
conditions becomes rather complicated as it re- 
quires taking into consideration several aspects 
arising from a particle leaving the computational 
domain from one periodic boundary has to reen- 
ter the domain back from a corresponding sec- 
ond periodic boundary with appropriate condi- 
tions. The periodic boundary conditions are im- 
plemented with the help of some appropriately 
defined transformation matrices. 

• The symmetric boundary condition is imple- 
mented in such a way to satisfy the criterion that 
for every particle crosssing the symmetry line, a 
similar one re-enters the domain in a direction 
given by the reflection off of the symmetry line. 

• When the particles try to move out of the exit 
boundary, they are removed out of the compu- 
tation. 

13 DETAILS OF THE COUPLING 
BETWEEN LSPRAY-II AND THE OTHER 

SOLVERS 

• The spray code is designed to be a stand-alone 
module such that it could easily be coupled with 
any other unstructured-grid CFD solver and the 
same holds true for EUPDF. (However, some 
grid-related parameters such as area vectors, grid 
connectivity parameters, etc. need to be sup- 
plied separately). 

• The spray solver needs information on the gas 
velocity and scalar fields from the other solvers 
and, then, it in turn supplies the liquid-phase 
source terms. 
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SPAWN MULTIPLE PEs FOR PARALLEL COMPUTING 


READ PARAMETERS & GEOMETRIC DATA 


[ CALL SPRAY PDF READ PARAMETERS 1 


INITIALIZE^ J_ 

H IMtTaT|ZE^DF — k/ARJABIUES 


CALL PDF INT RERUN 


RERUN 


READ PDF RESTART FILES 


INITIALIZE 


CALL SPRAY INT RERUN 




RERUN 


1 INITIALIZE SPRAY VARIABLES I [READ SPRAY RESTART FILES "1 

L . i I ' 


INITIALIZE 



RERUN 



OUTPUT 

I 

STOP 

Fig. 8 The overall flow structure of the combined CFD, LSPRAY-II, and EUPDF solvers. 
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Table 8. ncciiquid-solver.in file. 

Input file content 

comments 

heading 

title of controlling parameters 

ldread ,ispray _mod , 
(when_start_spray (n_i) , 
nd=l,no_of_injectors) 

If ldread = .true., restarts the calculation from the data stored from 
the previous computation. Otherwise initiates a new spray computation. 

ispray_mod controls the calls to the spray solver. The spray solver is called 
once at every other CFD iteration as given by the number, ispray_mod. 

when_start_spray() represents the CFD iteration where the computations for 
the n_ith injector are initiated. 

heading 

title of controlling parameters 

dtml, dtgl, & 
dtil 

dtml = time step for advancing the liquid phase equations. 

dtgl = global time step. Whenever the spray solver is called, it advances 
the spray computations over a period of dtgl before returning control 
over to the calling routine. To be more precise, it advances the 
liquid phase equations over a number of time steps as determined by 
(dtgl/dtml). 

dtil = injection time step. It is where a new group of droplets are 
introduced into the computation. 


• The PDF solver needs information on the mean 
gas velocity, turbulent diffusivity and frequency 
from the CFD solver and the liquid-phase source 
terms from the spray solver, and then it in turn 
provides the solution for the scalar (species and 
energy) fields to the flow and spray solvers. 

• It should also be noted that both the PDF and 
spray solvers are called once at every other spec- 
ified number of CFD iterations. 

• All of the three solvers (LSPRAY-II, EUPDF, 
and CFD) are advanced sequentially in an iter- 
ative manner until a converged solution is ob- 
tained. 

• All three codes (EUPDF, CFD, and LSPRAY- 
II) were coupled and parallelized in such a way 
to achieve maximum efficiency. 

The coupling issues could be better understood 
through the use of a flow chart shown in Fig. 8. 
It shows the overall flow structure of the combined 
CFD, LSPRAY-II, and EUPDF modules. Both the 
PDF and spray codes are loosely coupled with the 
CFD code. The spray code is designed in such a 
way that only a minimal amount of effort is needed 


for its coupling with the flow and PDF solvers. The 
present version of the spray module relies entirely on 
the use of Fortran common blocks for its informa- 
tion exchange with other modules. Even this reliance 
should entail only few changes to be made within the 
spray code for its linkage with different solvers. The 
PDF code is also structured along similar coupling 
principles. 

The flow chart of Fig. 8 contains several blocks 
- some shown in black and/or solid lines and the oth- 
ers in color and/or dashed lines. The ones in solid 
blocks represent the flow chart that is typical of a 
CFD solver. The ones in dashed blocks represent the 
additions arising from the coupling of the spray and 
PDF solvers. 

The coupling starts with the calling of the sub- 
routine - spray _pdf_read_parameters, which then 
reads the spray control parameters from the input file, 
nccdiquidjsolver.in of Table 8. This table provides 
a detailed description of the following input file vari- 
ables - ldread, isprayjmod, (when_start_spray(ni), 
nJ=l,no_of_injectors), dtml, dtgl, & dtil. The 
coupling is then followed by the calling of the 
pdf _int_rerun subroutine. It initializes PDF com- 
putations and, also, it may restart the PDF compu- 
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Fig. 9a An illustration of the parallelization 
strategy employed in the gas flow comp- 
utations . 
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Fig. 9b An illustration of the parallelization 
strategy employed in the spray computa- 
tions . 
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tations if needed from the data stored from a previ- 
ous iteration. Similarly, we call spr ay _int .rerun for 
the spray computations. It is noteworthy that the 
spray computations can be restarted by reading the 
data from the restart files - nccJiquid-params.out 
& nccJiquidjresults Ab. The restart capability is in- 
voked by setting the logical variable, Idread, of the 
input file, nccJiquidsolver.in of Table 8, to be true. 
Otherwise, the spray computations are initialized to 
start from the beginning. Then, the coupling pro- 
ceeds with the calling of the following subroutines: 
dclr for integrating the spray calculations and eupdf 
for the Monte Carlo PDF. The input variable, is- 
pray_mod of Table 8, controls the calls to the spray 
integrating routine. The spray solver is called once at 
every other number of CFD iterations as specified by 
ispray_mod. And the first call to the spray solver is 
controlled by the input variable, when_start_spray() 
of Table 8, which represents the starting CFD itera- 
tion number from which the spray computations are 
initiated. Finally, the coupling ends with the call- 
ing of a subroutine, spray_pdf .output, which will 
create a set of new restart files. 

14 DETAILS OF THE FORTRAN 
SUBROUTINES & FUNCTIONS 

Table 9 provides a list of all the Fortran sub- 
routines developed as a part of the spray module. 
This table also provides information on all the For- 
tran functions. It also describes the purpose of all 
the individual subroutines and functions. 

15 DETAILS OF THE PARALLELIZATION 

There are several issues associated with the par- 
allelization of both the spray & PDF computations. 
The goal of the parallel implementation is to extract 
maximum parallelism so as to minimize the execu- 
tion time for a given application on a specified num- 
ber of processors [37] . Several types of overhead costs 
are associated with parallel implementation which in- 
clude data dependency, communication, load imbal- 
ance, arithmetic, and memory overheads. The term 
arithmetic overhead is the extra arithmetic opera- 
tions required by the parallel implementation. Mem- 
ory overhead refers to the extra memory needed. Ex- 
cessive memory overhead reduces the size of a prob- 
lem that can be run on a given system and the 
other overheads result in performance degradation 
[37]. Any given application usually consists of sev- 
eral different phases that must be performed in cer- 
tain sequential order. The degree of parallelism and 


data dependencies associated with each of the sub- 
tasks can vary widely [37]. The goal is to achieve 
maximum efficiency with a reasonable programming 
effort [37]. 

In our earlier work, we discussed the parallel 
implementation of a spray algorithm developed for 
the structured grid calculations on a Cray T3D [10]. 
These computations were performed in conjunction 
the Monte Carlo PDF method. The parallel algo- 
rithm made use of the shared memory constructs ex- 
clusive to Cray MPP (Massively Parallel Processing) 
Fortran and the computations showed a reasonable 
degree of parallel performance when they were per- 
formed on a NASA LeRC Cray T3D with the number 
of processors ranging between 8 to 32 [10]. Later on, 
the extension of this method to unstructured grids 
and parallel computing written in Fortran 77 with 
PVM or MPI calls was reported in [11-15]. The lat- 
est version in Fortran 77/90 offers greater computer 
platform independence. In this section, we only high- 
light some important aspects of parallelization from 
Refs. [10-15]. 

• Both the EUPDF and CFD modules are well 
suited for parallel implementation. For the gas- 
phase computations, the domain of computation 
is simply divided into n-Parts of nearly equal 
size and each part is solved by a different pro- 
cessor. Fig. 9a illustrates a simple example of 
the domain decomposition strategy adopted for 
the gas-phase computations where the total do- 
main is simply divided equally amongst the avail- 
able computer processing elements (PEs). In 
this case, we assumed the number of available 
PEs to be equal to four. 

• But the spray computations are more difficult to 
parallelize for the reasons summarized below: 

(1) Non-uniform nature of spray distribution: 
Most of the particles are usually confined to a 
small region where the atomizer is located. 

(2) Dynamic nature of Lagrangian particles: 
Particles keep moving between different sub- 
divided domains of an Eulerian grid (grid used 
in the CFD computation) which are assigned to 
different processors. While some new particles 
are introduced at the time of fuel injection, some 
others are taken out of computation. 

Conceptually, there are several ways to paral- 
lelize the spray computations, we, however, devel- 
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Table 9. Description of LSPRAY-II Fortran subroutines & functions. 

Function 

Purpose of the Function 

blasiu(x) 

This function returns a solution for the function, f{Bk ), 
of Eq. (21) for use in computing the droplet regression rate. 

peng_rob_cmpr 
(chem_model,ycomp, 
tij, density, 
number _of_species) 

It provides a value for the pressure based on the solution 
of the Peng- Robinson EOS. 

Subroutine 

Purpose of the Subroutines 

chaslv 

This routine has two main functions: 

(1) It integrates the liquid phase equations. 

(2) It removes the particles that need to be taken 
out of the computation. 

dclr 

This routine is called once at every other CFD iteration as 
specified by ispray_mod. It is primarily a controlling 
routine for spray computations. It is only called in 
conjunction with steady _state_model = .false. 

This routine has the following functions: 

(1) It initializes the source terms to zero. 

(2) Checks to see if new particles need to be introduced. 

(3) Advances liquid phase equations over an allowable or 
pre-specified time step, dtml, w r ith calls to the following 
routines: 

intplal - Interpolates the gas phase properties at the particle 
location. 

chaslv - Advances liquid phase equations. 

intpla - Identifies computational cells and PEs associated 
with particles. 

sprips - Evaluates the liquid phase source term contributions 
of the CFD and PDF equations. 
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Table 9. Description of LSPRAY-II Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 


(4) Continues with steps (2) and (3) until the computations 
are completed over one single global time step, dtgl. 

(5) Returns control over to other solvers, e.g. CFD and EUPDF, 
and supply them with the source terms averaged over dtgl. 

dclr .steady 

This routine is similar to dclr but used only in conjunction 
with steady _state_model = .true. 

This routine has the following functions: 

(1) It initializes the source terms to zero. 

(2) Introduces a new set of spray particles. 

(3) Advances liquid phase equations over an allowable or 
pre- specified time step, dtml, with calls to the 
following routines: 

intplal - Interpolates gas phase properties at the 
particle locations. 

chaslv - Advances liquid phase equations. 

intpla - Identifies computational cells and PEs associated 
with particles. 

sprips - Evaluates the liquid phase source term contributions 
of the CFD and PDF equations. 

(4) It repeats step (3) until there are no more spray particles 
left to be integrated. (NOTE: The particles are taken out of 
computation when they reach either a certain size of negligible 
proportion or exit out of the computational domain.) 

(5) Returns control over to other solvers, e.g. CFD and EUPDF, 
and also supply them with source terms. 
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Table 9. Description of LSPRAY-II Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 

dropdis(rhol, 

flowdum,sr, 

fld,smd,nofg) 

This routine computes droplet distribution from the 
following correlation (also, appears as Eq. 53): 

dn/n = a((D/D 32 ) al r)exp(-b((D/D 32 ) bet ))dD/D 32 
where a, b, alp, and bet are constants. 

find_transport_ds ( 
element, ijle,ymgf) 

It computes the following properties of the gas mixture 
at the droplet interface by making use of the one-third 
rule of Eq. (46): molecular viscosity, gas density, 
and diffusion coefficient. 

find_xyzface 

(i) 

This routine computes x, y, and z locations of all the face 
centers of an element, i. This information is used in the 
particle search algorithm. 

getiiq_tat_properties 

(tmp,pmp,j) 

Computes the following variable properties of a liquid 
mixture: density, specific heat, molecular viscosity, 
gas density, and diffusion coefficient. 

hp_transport 

(tij,pij,x_work, 

its,avisc,acond) 

It calculates the high pressure correction for the gas-phase 
transport properties: viscosity, conductivity, and diffusion. 

hp_diff_table 

(ct,pt,dabcvi) 

This subroutine provides a value for {Dab P)/(D^b P) + 
based on the Takahashi high-pressure correlation for the 
calculation of the diffusion coefficients in a gas mixture. 

For further details, see Reid et. al. [16], “The properties of 
Gases and Liquids.” The correlation data is read from 
the file - hp_diff_table.dat. 

intpla 

This routine performs three main functions: 

(1) Particle Tracking - It identifies the computational cell 
in which a particle is located. In parallel computing, 

it also means identifying the corresponding processor of the 
computational domain in which a particle is located. 

(2) It implements appropriate boundary conditions. 

(3) It reassigns the particles between different PEs based on 
the parallel strategy employed. 
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Table 9. Description of LSPRAY-II Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 

intplal 

This routine interpolates the gas-phase properties at the 
particle location. In the present case, a simple first-order 
interpolation is employed. 

mimd .spray 

For computational elements whose neighboring cells are 
assigned to a different processor, it initializes arrays, 
ipr_frud() and ile_frdd(), for storing information on the 
processor and element ID numbers of the neighboring cells. 
This information is needed in order to process the particle 
movement between the domains of the neighboring processors. 

mimd_spray_recv 

(i-recvfrom) 

This subroutine is called by mimd_spray in order to 
gather the relevant information from the neighboring 
processors. 

mimd_spray_send 

(i-sendto) 

This subroutine is also called by mimd-spray in order 
to send some relevant information to the neighboring 
processors. 

par_loc(xparz, 

yparz,zparz, 

ipare,iparp) 

Given the x,y,z coordinates of a particle location, the 
algorithm identifies the corresponding computational cell 
in which a particle is located. It also identifies the 
corresponding processor of the computational grid in which a 
particle is located. This information is useful when new 
droplet groups are introduced at the time of injection. 

peng_rob 

(chem_model,ycomp, 
tij pressure, 
number_of_species) 

It provides a solution based on the Peng-Robinson EOS which 
is applicable for real gases at high pressures. It provides a 
value for the density of a gas mixture as the output 
variable. 

peng_rob_gen 

(xsp,ppr,tpr, 

its,rho0) 

It takes the following as the input variables: xsp() - the mass 
fraction of the species, ppr - pressure, tpr - temperature, and 
its - the number of species. And it returns as the output, rhoO - 
density based on the solution of the Peng-Robinson EOS. 
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Table 9. Description of LSPRAY-II Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 

peng_rob_sc 

(ppr,tpr,its) 

It provides a solution for Z* c and pi r for ith species based 
on the solution of the Peng- Robinson EOS for a given value 
of Pi r and t i r . It also requires the specification of the 
total number of species, its. 

pr_eos_read_files 

This routine reads the following three input files: 

For ncc_liquid_physical_properties.dat, it reads the 
following input variables of Table 1: 

tboil: normal boiling point, K 

tcrit: critical temperature, K 

pcrit: critical pressure, atm 

rholn: liquid density (at 1 bar, 273.15 K), kg/m 3 

elhin: heat of vaporization at normal boiling point, KJ/kg 

vole: critical volume, cm 3 /g-mole 

volln: molar volume at normal pressure, cm 3 /g-mole 

aom: Pitzer’s acentric factor 

sigma: characteristic diameter of a molecule, angstroms A 0 
ep_ka: epsilon /kappa: K 

For ncc_liquid_binary_parameters.dat, it reads the values 
of the binary parameters, K of Table 2. 

For nccJiquid_transport.dat, it reads the polynomial 
polynomial constants of Tables 3-5 used in the evaluation 
of C p , thermal conductivity, and viscosity. 
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Table 9. Description of LSPRAY-II Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 

prnspr 

It writes standard spray output data. 

root(idf,a,b,c,zv,zl) 

It provides a value for the compressibility factor, Z, from 
the solution of Eq. (6). One of the three roots of Eq. (6) 
yields Z v and the other Z/. 

spr ay _int -rerun 

This routine initializes the spray solver based on the input 
read from various input files containing different spray 
parameters. In addition to initializing some of the variables, 
it also restarts the computations from the stored data of a 
previous computation if it is a rerun. 

spray _plot- 
_output 

It writes some useful plotting data when used in conjunction 
with the steady_spray .model = .false. 

sprips 

This routine computes the liquid-phase source terms of Eqs. 
(l)-(4) for use in both CFD and Monte Carlo PDF solvers. 

smlc(i) = liquid-phase contribution of Eq. (1) of Section 
4. 

smlmx(i), smlmy(i), smlmz(i) = liquid-phase contribution 
of Eq. (3) of Section 4. 

smle(i) = liquid-phase contribution of Eq. (4) of Section 4. 

sy(il,iu,bb, 

dd,aa,cc) 

It is a tri-diagonal matrix solver. It is used in the solution 
of both Eqs. (27) & (37). 

uvw_par(swlr_angle(nx) , 

angle_work,nx,ny,v_inj , 

nmip,t_rotation, 

cone,n_cone_rays, 

cone_rotation,uloc, 

vloc,wloc,nJ) 

It computes the particle injection velocity for different 
cone configurations of Figs. 2 to 4. 


Table 10. CPU time (sec) per cycle versus number of PEs. 



Number of processors 

Solver 

Characteristic 

2 

5 

10 

CFD 

5 steps/cycle 

2.50 

1.25 

0.75 

EUPDF 

1 step/ cycle 

6.5 

2.9 

1.9 

LSPRAY-II 

100 steps/cycle 

1.70 

0.64 

0.53 
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Table 11. CPU Time (Sec) Per Cycle Versus Number of PEs. 


Number of Processors 

Solver 

Characteristic 

1 

2 

4 

8 

16 

Spray 

100 steps/cycle 

6.83 

5.29 




Max. Spray Particles in a PE 

2695 

2097 

1165 

623 

312 

Min. Spray Particles in a PE 

2695 

598 

118 

14 

0 


oped and tested two different domain decomposition 
strategies [10-14]. 

• Strategy I: 

The Lagrangian particles were assigned fairly 
uniform amongst the available processors but the 
calculations associated with the particle track- 
ing, the interpolation of the gas-phase proper- 
ties, and the source-term evaluation were com- 
puted on the processor of the computational grid 
in which a particle is located. 

This strategy leads to an uniform loading during 
integration but leads to excessive message pass- 
ing. 

• Strategy II: 

The Lagrangian particles were assigned to the 
processor of the computational grid where the 
particle is located. Fig. 9b illustrates a sim- 
ple example of the domain decomposition strat- 
egy adopted for the liquid-phase computations 
where the corresponding gas-flow computational 
domain is divided into equal parts between the 
four available PEs. In this strategy, the La- 
grangian particles are assigned to the processor 
of the computational grid where a particle is lo- 
cated. 

This strategy lead to a non-uniform loading dur- 
ing integration but leads to less message passing. 

Our experience has shown that Strategy II 
seems to work well on different computer platforms: 
both massively parallel computers as well as hetero- 
geneous cluster of workstations. So in the present ver- 
sion of the code, we have opted to implement Strategy 
II over Strategy I. 


16 DETAILS OF PARALLEL 
PERFORMANCE 

The details of the combined parallel perfor- 
mance of the CFD, EUPDF, and LSPRAY-II codes 
involving several different cases can be found in Refs. 
[10-15]. Here, we only summarize briefly the parallel- 
performance results for two different cases. One is 
a 3D test case and more details on this case can be 
found in the reference [13]. For this case, the calcu- 
lations were performed on a computational grid com- 
prising of 8430 tetrahedral elements and 100 Monte 
Carlo PDF particles per cell. The computations were 
performed on one of the NASA Ames Research Cen- 
ter’s parallel computer platforms called Turing which 
is a SGI Origin work-station with 24 PEs (Processor 
Elements). Table 10 summarizes the CPU times per 
cycle taken by the EUPDF, LSPRAY-II, and CFD 
solvers vs the number of PEs. Both the CFD and 
PDF solvers show good parallel performance with an 
increase in the number of processors but for the spray 
solver it shows reasonable parallel performance. 

Next, we would like to summarize the results 
from [13, 38] showing only the results of spray compu- 
tations. The results are summarized in Table 11. The 
computations were performed on Turing at NASA 
ARC (Ames Research Center) - it is a SGI Origin 
work-station with a maximum of 24 PEs (Processor 
Elements). The computations made use of an un- 
structured grid with a mesh size of 3600 tetrahedral 
elements. And it made use of about 2,695 Lagrangian 
particles for the spray computations and one hundred 
Monte Carlo particles per element for the PDF com- 
putations. In a given cycle of one global time step, 
dt g i , the spray equations were advanced over one hun- 
dred time steps as given by dt g i/dt m i = 100. The 
first row of Table 11 summarizes the CPU times per 
PE per cycle taken by the spray solver vs the num- 
ber of PEs. As expected, the CPU time goes down 
with an increase in the number of processors. In an 
ideal case, one would expect an inverse reduction in 
epu time with an increase in the number of proces- 
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sors. Here, we don’t get such an ideal performance 
because of the resulting non-uniform distribution of 
spray particles, between various participating proces- 
sors, from the implementation of Stratergy II. To get 
an idea of the spray particle distribution, we have tab- 
ulated the maximum and minimum number of parti- 
cles found between various processors. When we go 
from 1 to 2 PEs, 2097 particles are assigned to one 
and the rest to the second. With four they are dis- 
tributed between 1165 and 118, with eight between 
623 and 14, and with sixteen from 312 to 0. The 
results clearly show that the reduction in the CPU 
time varies almost linearly with the reduction in the 
number of maximum particles. 

17 A SUMMARY OF SOME RECENT 
VALIDATION CASES INVOLVING BOTH 
REACTING AND NON-REACTING 
SPRAY COMPUTATIONS 

A total of the following five cases were validated: 

1. A reacting methanol spray with no-swirl. 

2. A non-reacting methanol spray with no-swirl. 

3. A confined swirl-stabilized n-heptane reacting 
spray. 

4. An unconfined swirl-stabilized n-heptane react- 
ing spray. 

5. A confined swirl-stabilized kerosene reacting 
spray. 

The experimental data for the first two cases 
was provided by McDonell & Samuelsen from the 
University of California at Irvine [38]. Both the cases 
are without swirl; one is a reacting case and the other 
is non-reacting. The data for the third and fourth 
cases was provided by Bulzan from the NASA Glenn 
Research Center [39-40]. Both the cases are swirl- 
stabilized reacting cases, one is an unconfined flame 
and the other is confined. The data for the last case 
was provided by El Banhawy & Whitelaw from Im- 
perial College [4]. It is a confined swirl-stabilized 
kerosene spray flame. The first three cases made 
use of unstructured grids and the last two structured 
grids. 

Here, we would like to provide a brief summary 
of the validation cases but a detailed presentation of 
the results and discussion can be found elsewhere in 
the papers [10-13]. The comparisons involved both 


gas and drop velocities, drop size distributions, drop 
spreading rates, and gas temperatures. The results 
were in reasonable agreement with the available ex- 
perimental data. The comparisons also involved the 
results obtained from the use of the Monte Carlo PDF 
method as well as those obtained from a conventional 
CFD solution without the Monte Carlo PDF method. 
For the first case of McDonell & Samuelsen’s reacting 
spray flame, the detailed comparisons clearly high- 
lighted the importance of chemistry /turbulence in- 
teractions in the modeling of reacting sprays [13]. 
The results from the PDF and non-PDF methods 
were found to be markedly different with the PDF 
solution providing a better approximation to the re- 
ported experimental data. The PDF solution showed 
that most of the combustion occured in a predom- 
inantly diffusion-type of flame environment and the 
rest occuring in a predominantly premixed-type of 
flame environment. However, the non-PDF predic- 
tions showed incorrectly that most of the combustion 
occured in a predominantly vaporization-controlled 
regime. The Monte Carlo temperature distribution 
showed that the functional form of the PDF for the 
temperature fluctuations varied substantially from 
point to point. The results brought to the fore some 
of the deficiencies associated with the use of assumed- 
shape PDF methods in spray computations. 

18 CONCLUDING REMARKS 

• This manual provides a complete description of 
LSPRAY-II - a Lagrangian spray solver devel- 
oped for application with parallel computing and 
unstructured grids. 

• It facilitates the calculation of the multi- 
component liquid sprays with variable properties 
valid over a wide range of low pressure condi- 
tions. 

• It provides the user with a basic understanding 
of the the spray formulation and the LSPRAY-II 
code structure, and complete details on how to 
couple the spray code to any other flow code. 

• The basic structure adopted for the grid repre- 
sentation and parallelization for the gas side of 
the flow computations follows the guidelines es- 
tablished for NCC. 

• Also, we have extended the joint scalar Monte 
Carlo PDF method to two-phase flows and, 
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thereby, demonstrating the importance of chem- 
istry/turbulence interactions in the modeling of 
reacting sprays. 

• Based on the validation studies involving several 
confined and unconfined spray flames, the results 
were found to be encouraging in terms of their 
ability to capture the overall structure of a spray 
flame. 

• The source code of LSPRAY-II will be available 
with NCC as a complete package. 
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